VELOCITY CORRELATIONS IN 

DRIVEN TWO-DIMENSIONAL GRANULAR MEDIA^ 



C. BIZON, M. D. SHATTUCK, 

J. B. SWIFT AND HARRY L. SWINNEY 

Center for Nonlinear Dynamics and Dept. of Physics, 
University of Texas, Austin, TX 78712, USA 



Abstract. Simulations of volumetrically forced granular media in two di- 
mensions produce states with nearly homogeneous density. In these states, 
long-range velocity correlations with a characteristic vortex structure de- 
velop; given sufficient time, the correlations fill the entire simulated area. 
These velocity correlations reduce the rate and violence of collisions, so 
that pressure is smaller for driven inelastic particles than for undriven elas- 
tic particles in the same thermodynamic state. As the simulation box size 
increases, the effects of velocity correlations on the pressure are enhanced 
rather than reduced. 



1. Introduction 

In rapid flows of granular media, the mean time between collisions of grains 
is much longer than the duration of a collision [1]; for such flows, the ma- 
chinery of kinetic theory is expected to apply. Continuum equations [2, 3] 
analogous to the Navier-Stokes equations can be produced, allowing quanti- 
tative analysis of flows. The simplest and most common formulations incor- 
porate Boltzmann's assumption of molecular chaos: that particle velocities 
are uncor related. 

While this assumption works well for low-density molecular gases, gran- 
ular gases may not abide such a restriction because collisions between 
grains are inelastic. Inelastic collisions reduce relative velocities, so that 
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post-collisional velocities are more parallel than pre-collisional velocities. 
Repeated inelastic collisions can lead to strong, long-range velocity correla- 
tions, which standard kinetic theory does not include. We will use molecular 
dynamics simulations to produce steady state granular gases and study the 
velocity correlations that develop. 

The importance and intrinsic interest of velocity correlations in gran- 
ular flows have been noted by a number of researchers. Two-dimensional 
simulations of an initially homogeneous distribution of inelastic disks with- 
out velocity correlations show that as time progresses, velocity correlations 
build in both strength and range [4] . These simulations are limited in time, 
however, because the homogeneous state is unstable to density fluctua- 
tions, and rapidly becomes inhomogeneous. Nevertheless, these simulations 
clearly displayed a characteristic vortex structure of the correlations. Based 
upon similar considerations, ring kinetic theory, which accounts for veloc- 
ity correlations, has been applied to the cooling state [5]. One-dimensional 
simulations of stochastically forced point particles also show velocity cor- 
relations [6]. 

We apply stochastic forcing [6, 7] to two-dimensional event-driven simu- 
lations of inelastic disks. The forcing overcomes the tendency of the granular 
material to form density clusters, and approximately homogeneous steady 
states form. In an earlier study of these states [8] , we found strong velocity 
correlations that extended throughout the entire simulation area. In the 
present work, we discuss the simulation method, show that the velocity 
correlations are essentially independent of the simulated area, and describe 
the vortex structure of the correlations. 

2. Simulations of Driven Granular Gases 

We treat collisions between molecules as instantaneous and binary. The 
collisions between grains conserve momentum but dissipate energy. Between 
collisions, particles travel along straight lines if unaccelerated, or along 
parabolas if accelerated. This model allows efficient simulation of collections 
of particles using event-driven molecular dynamics [9, 10]. 

When particles collide, the component of the relative particle velocity 
along the line joining particle centers, Vn, is reversed, and reduced by a 
factor e, the coefficient of restitution, which can take values between 1 
for elastic particles and for completely inelastic particles. We allow e to 
depend on Vn through 




,Vn<Vo 
,Vn> Vo 



(1) 
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where B = {I - e)ivo)~f^, P = 3/4 and e is a constant, chosen to be 0.7. 
These parameters give quantitative agreement to experiments on patterns 
in vertically oscillated granular media [11, 12]. The variation in e has the ef- 
fect of removing inelastic collapse [13], which is a singularity in the inelastic 
hard sphere model that produces an infinite number of collisions within a 
finite time [14, 15]. In general, colliding particles also exert frictional forces 
on one another; for this paper, we assume that the coefficient of friction is 
zero, so that we are studying only the effects of inelasticity. 

Because of inelasticity, the energy of an unforced collection of grains 
inevitably decreases. To achieve steady states, then, we must force the 
granular material. Methods that force through boundaries, such as shak- 
ing, invariably produce strong inhomogeneities in the system; to achieve 
near-homogeneity, we force volumetrically, assuming the particles to be in 
contact with a white-noise heat bath [7]. Whenever two particles collide, the 
velocities of two other randomly selected particles are changed by amounts 
|5v|fj, where the magnitude of the kicks, |(5v|, are always the same, but 
the direction vectors, fj arc randomly chosen for each kicked particle. In 
addition to the white noise heat bath, we perform a lesser number of runs 
with two other heat baths. To model the motions of pucks on an air ta- 
ble [16, 17], we can allow particles to accelerate randomly from collision 
to collision. Finally, we model the effects of a strong heat bath, which we 
denote the Boltzmann bath, by completely obliterating the velocities of 
randomly chosen particles, and giving new velocities based on a Boltzmann 
distribution. The details of all three forcing methods may be found in [8]. 

We perform simulations of N disks of diameter a moving in a two- 
dimension square of side length L, which varies from 52.6a to 420. 8c7. The 
simulation box is periodic in both directions. The solid fraction, defined 

2 

as Njj^, is 0.5 for all runs. Because of the variation of e with relative 
normal velocity, the velocity scale vq enters; we use vo to nondimensionalize 
velocities, and Vq to nondimensionalize the granular temperature T. For T 
much larger than one, most particle collisions will occur with the high- 
velocity value of e, 0.7; for lower T, a range of e will occur. 

3. Dependence of Correlations upon Simulation Area 

We denote two particles 1 and 2, and k the a unit vector pointing from 
the center of 1 to the center of 2. The velocity of 1 then has a components 
parallel to, and perpendicular to, v^, k, as does particle 2. We define 
two correlation functions 




(2) 




r-i 



(3) 
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where the sums are over the particles such that the distance between the 
two particles is within Sr of r. For uncorrelated particle velocities, (^1^2) 
and {V1V2) will both give zero. 

In the smallest simulation area, L = 52. 60", correlations extend the full 
length of the computational cell. Cell filling structures may be divided into 
two cases: structures with a natural length that is larger than the box 
in which they exist and structures that will always grow to fill any finite 
box. To differentiate between the former and the latter, we performed four 
simulations with white noise forcing, quadrupling the area at each step, 
while holding the solid fraction fixed at 0.5. The granular temperature T is 
approximately 30, but varies between 28 in the smallest box and 32 in the 
largest. This variation in temperature is not important; for T >> 1, the 
coefficient of restitution is independent of collision velocity. In this limit, 
the role of the temperature is simply to set the velocity scale. The velocity 
correlation functions are shown in Fig. 1. Even in the largest simulation, 
composed of 112768 particles, the correlations fill the box. However, the 
correlation functions for the largest simulation are somewhat different from 
the smaller ones. This is probably due to poorer statistics; in terms of 
collisions per particle, this run lasted only one-half as long as the next 
largest. 

Because velocity correlations are positive for small separations, particles 
collide less frequently and with less relative velocity than elastic particles 
at the same density, for which velocity correlations are much smaller. As a 
result, less momentum will be transferred through inelastic collisions than 
through elastic collisions, and the pressure, P, will decrease. 

Assuming that velocity correlations do not exist, the equation of state 
for dense granular gases is given by [3] 

P={A/TTa^)iyT{l + {l + e)G{iy))., (4) 

The first term on the right hand side, (4/7ro"^)j^T, accounts for momentum 
transfer due to particle streaming without collisions, while the second term, 
{4/Tra'^)uT{l + e)G{i'), accounts for the momentum transfer due to particle 
collisions [18]. In the absence of velocity correlations, G(z/) is defined as 
vgiy^a), where 5(1^,0") is the radial distribution function for the particles, 
evaluated at zero particle separation. Calculation of P from simulation, 
via measurement of the virial [19], becomes a measurement of G(z^), which 
describes the collisional momentum transport. If velocity correlations exist, 
G{i') will be reduced, since less momentum will be transported collisionally. 

Figure 1 shows that the short range velocity correlations depend on the 
size of the box; therefore, G{u) should also depend on L. Figure 2 displays 
G{i') as a function of L for these four runs. Over about one decade, G{i') 
scales with logL. Clearly this scaling can not continue indefinitely, since 
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Figure 1, Velocity correlations as a function of particle separation at = 0.5 and T « 30, 
for four different box sizes. + : L = 52a, A : L = 105(7, * : L = 211a, a : L = 421c7. 



unphysical negative values of G{i') would result. Note also, that increasing 
the box size actually leads to values of G{u) farther from the values for 
uncorrelated velocities. 

This unusual result, that the importance of velocity correlations in- 
creases with increasing computational area, can also be deduced from the 
distribution of collision velocities. Figure 3 exhibits these distributions for 
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0.0 0.5 1.0 1.5 2.0 2.5 3.0 
log.(L/U) 



Figure 2. G as a function of L for the runs shown in Fig. 1. Lo = 52a denotes the length 
of the smallest box. The dotted line is a fit to all four points: G{iy) = 1.3—0.04 log2(i/io), 
The solid line as a fit to the three largest L values: G(i') = 1.295 — 0.038 log2(ij/ijo)- 
Note that the log is base 2. 

the runs displayed in Figures 1 and 2. As the computational area increases, 
so too does the deviation from the distribution predicted for particles cho- 
sen without correlation from a Boltzmann distribution, plotted as a solid 
curve. 

4. Vortex Structure 

Inelasticity breeds velocity correlations; reduction of relative velocity in 
collisions leads to particles moving more alike after collisions than before. 
On average, then, particles will be surrounded by particles that are moving 
along with them. The structure of the velocity correlations can be eluci- 
dated by calculating this average flow around each particle. 

For a single particle i, we can calculate the flow around it by translating 
it to the origin, and rotating so that its velocity lies along the positive x 
axis. If v(a;, y) is the velocity field defined by the particles, then the flow 
around particle i is given by 

Ui = Re{i)^ix -Xi,y- yi), (5) 

where 6{i) is the angle between the i-th particle velocity, Vj , and the 
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Figure 3. Probability distribution of collision velocities Vc = |vi — V2I, for the data in 
Figs. 1 and 2. + : L = 52cr, A : L = lOStr, * : L = 211a, □ : L = 421.7. The solid curve 
is P{vc/VT) = (l/2\/7rT3)t;2e-"c/4T^ ^iij^h holds for elastic particles. 

positive X axis, {xi,yi) is the position of the i-th particle, and Rg is the 
operator that rotates vectors clockwise through angle 9. The average flow 
around particles, then, is 

TV 

u = ^u,/iV. (6) 

1=1 

Finally, u is averaged over about 100 frames to reduce noise. 

Figure 4 displays vector fields of the average flow around particles, u, 
for the three types of forcing, as well as for unforced elastic particles, all 
at u = 0.5 and T = 1.05. In each case, the vector at the origin, which 
measures only the average particle speed, has been suppressed, and the 
longest remaining vector in each field has been scaled to unit length. In both 
the white noise and accelerated forcings, the average flow near the origin 
is along the positive x axis, i.e., with the direction of the central particle's 
motion. The Boltzmann bath shows some indications of this effect close to 
the origin, but the correlations are destroyed by the strongly thermalizing 
forcing before they can propagate to larger length scale. For the elastic 
particles, there is no discernible flow, only noise. 

Close to any particle, surrounding particles move along with it. Far- 
ther away, the correlations decay and cannot be seen on Fig 4, so the 
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Figure 4 ■ The average velocity fields around a particle centered in each cell and moving 
to the right, u, for elastic particles and for inelastic particles forced in three different 
ways (cf section 2). Each vector field is scaled separately so that its longest vector has 
length one. Compared to the (suppressed) central vector, these lengths are: White noise, 
0.2; Accelerated, 0.27; Boltzmann, 0.008; Elastic 0.008. The boxed regions in the white 
noise and elastic flows are shown in Fig. 5. 



boxed regions for the white noise forcing and for elastic particles are ex- 
panded in Fig. 5. While expansion of the velocity field for elastic particles 
produces still more noise, the inelastic flow field reveals a highly ordered 
vortex structure. Along the direction of the central particle's motion, the 
velocities slowly drop to zero, while perpendicular to the original particle's 
motion, the velocities drop to zero and increase in the negative direction; 
this flow makes clear the structure of the velocity correlation functions in 
Fig. 1. 
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Figure 5. A close-up on the boxed regions in Fig 4 reveals that for inelastic particles, 
large vortices form, one on each side of the particle. The longest vector in the velocity 
field for inelastic particles rcprosouts a velocity nine times larger than that represented 
by the longest vector for elastic particles. 
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This vortical flow is reminiscent of similar structures produced in simu- 
lations of elastic particles [20, 21] by Alder and Wainwright. In their simu- 
lations, they discovered diffusive behavior different from that predicted by 
kinetic theory. The diffusion constant may be written in terms of the slope 
of the exponentially decaying autocorrelation function. However, Alder and 
Wainwright found deviations from exponential decay, and traced the devi- 
ations to a vortical flow. If particles a and b are initially uncorrelated, an 
elastic collision will correlate each particle's post collision velocity with the 
other particle's pre-collision velocity; both particles now have a correlation 
with the original velocity of particle a. As particle b collides with other 
particles, they gain information about particle a's initial velocity. Several 
collision times later, this information has been transmitted to many parti- 
cles. 

There arc two main differences between the vortices in flows of elastic 
particles and those in flows of inelastic particles. Alder and Wainwright 
produced the flow field given by 

u{t')i = Re^i,t)^{x -Xi{t),y- yi{t),t'). (7) 

For t' = t, Eq. 5 is recovered; for clastic particles, no structure is apparent. It 
is only at later times, t' > t, that a vortex appears in u(t'). For the inelastic 
particles, however, structure is clear at t' = t. The second difference is the 
strength of the vortex. The strongest velocity in Alder and Wainwright 's 
vortex was about 2% of the original velocity, while for inelastic particles, 
the strongest velocity can be about 40% of the central velocity. 

The inelastic vortex is so strong that hints of it are visible even in a 
single snapshot of particles. Figure 6 shows such a snapshot, with particles 
colored black if they have positive horizontal velocity and white if they have 
a negative horizontal velocity. For elastic particles, the black and white are 
well mixed, but in the inelastic case larger scale structure can be glimpsed. 
Black particles are concentrated along the top and bottom of the image, 
and white particles are concentrated along the central region. 

5. Conclusion 

The correlations we have found are consistent with those of simulations 
on the homogeneous cooling state [4]. In those simulations, the range of 
velocity correlations grew until the onset of large scale density variations. 
The addition of forcing in our simulations suppresses the growth of density 
fluctuations, allowing the velocity correlations to continue to grow until 
they extend throughout the entire computational area. 

The results we obtain are not particularly sensitive to the exact form 
of the forcing. In both the white noise and accelerated forcing schemes. 
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Figure 6. Snapshots of simulations with white noise driving and with elastic particles; 
large coherent structures are visible for the dissipative system on the left. Particles with 
positive horizontal velocity are black, particles with negative horizontal velocity are white, 
(i^ = 0.5, T = 1.05) 



vortical correlation structures form. Only when the bath explicitly destroys 
correlations, as in the Boltzmann bath, do the results differ. 



6. Acknowledgments 

We thank J. T. Jenkins, M. H. Ernst, and T. P. C. van Noije for useful 
discussions. This work was supported by the Department of Energy Office 
of Basic Energy Sciences. 



References 

1. C. S. Campbell, Annu. Rev. Fluid Moch. 2, 57 (1990). 

2. C. K. K. Lun, S. B. Savage, D. J. Jeffrey, and N. Chepurniy, J. Fluid Mech. 140, 
223 (1983). 

3. J. T. Jenkins and M. W. Richman, Arch. Rat. Mech. Anal. 87, 355 (1985). 

4. J. A. G. Orza, R. Brito, T. P. C. van Noije, and M. H. Ernst, Int. J. Mod. Phys. C 

8, 953 (1998). 

5. T. P. C. van Noije, M. H. Ernst, and R. Brito, Phys. Rev. E 57, R4891 (1998). 

6. M. R. Swift, M. Boamfa, S. J. Cornell, and A. Maritan, Phys. Rev. Lett. 80, 4410 

(1998). 

7. D. R. M. Williams and F. C. MacKintosh, Phys. Rev. E 54, R9 (1996). 

8. C. Bizon, M. D. Shattuck, J. B. Swift, and H. L. Swinney, Submitted to Phys. Rev. 
E (1999). 

9. B. D. Lubachevsky, J. Comp. Phys. 94, 255 (1991). 

10. M. Marm, D. Risso, and P. Cordero, J. Comput. Phys. 109, 306 (1993). 

11. C. Bizon, M. D. Shattuck, J. B. Swift, W. D. McCormick, and H. L. Swinney, Phys. 
Rev. Lett. 80, 57 (1998). 



12 



12. J. R. do Bruyn, C. Bizon, M. D. Shattuck, D. Goldman, J. B. Swift, and H. L. 
Swinncy, Phys. Rev. Lett. 81, 1421 (1998). 

13. D. Goldman, M. D. Shattuck, C. Bizon, W. D. McCormick, J. B. Swift, and H. L. 
Swinncy, Phys. Rev. E 57, 4831 (1998). 

14. S. McNamara and W. R. Young, Phys. Fluids A 4, 496 (1992). 

15. S. McNamara and W. R. Young, Phys. Rev. E 50, R28 (1994). 

16. L. Oger, C. Annie, D. Bideau, R. Dai, and S. B. Savage, J. Stat. Phys. 82, 1047 
(1996). 

17. I. Ippolito, C. Annie, J. Lemaitre, L. Oger, and D. Bideau, Phys. Rev. E 52, 2072 

(1995). 

18. S. Chapman and T. G. Cowling, The Mathematical Theory of Non-uniform Gases 
(Cambridge University Press, London, 1970). 

19. D. C. Rapaport, The Art of Molecular Dynamics Simulation (Cambridge University 
Press, Cambridge, 1980). 

20. B. J. Alder and T. E. Wainwright, J. Phys. Soc. Japan 26, 267 (1968). 

21. B. J. Alder and T. E. Wainwright, Phys. Rev. A 1, 18 (1970). 



